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We enumerate and classify nearly all of the possible mechanically stable (MS) packings of 
bidipserse mixtures of frictionless disks in small sheared systems. We find that MS packings form 
continuous geometrical families, where each family is defined by its particular network of particle 
contacts. We also monitor the dynamics of MS packings along geometrical families by applying qua- 
sistatic simple shear strain at zero pressure. For small numbers of particles (N < 16), we find that 
the dynamics is deterministic and highly contracting. That is, if the system is initialized in a MS 
packing at a given shear strain, it will quickly lock into a periodic orbit at subsequent shear strain, 
and therefore sample only a very small fraction of the possible MS packings in steady state. In 
studies with N > 16, we observe an increase in the period and random splittings of the trajectories 
caused by bifurcations in configuration space. We argue that the ratio of the splitting and contrac- 
tion rates in large systems will determine the distribution of MS-packing geometrical families visited 
in steady-state. This work is part of our long-term research program to develop a master-equation 
formalism to describe macroscopic slowly driven granular systems in terms of collections of small 
subsystems. 



PACS numbers: 81.05.Rm, 83.80.Fg, 83.80.Iz 



I. INTRODUCTION 

Dry granular materials are collections of discrete, 
macroscopic particles that interact via dissipative and 
purely repulsive interactions, which are nonzero when 
particles are in contact and vanish otherwise. Granu- 
lar systems range from model systems composed of glass 
beads to pharmaceutical powders, to soils and geological 
materials. 

A distinguishing feature of granular materials is that 
they are athermal. Since individual grains are large, ther- 
mal energy at room temperature T is unable to displace 
individual grains. Thus, without external driving, gran- 
ular materials are static and remain trapped in a single 
mechanically stable (MS) grain packing with force and 
torque balance on each grain. In contrast, when external 
forces are applied to granular materials, these systems 
flow, which gives rise to grain rearrangements, fluctua- 
tions in physical quantities like shear stress and pressure, 
and the ability to explore configuration space. 

There are many driving mechanisms that generate 
dense flows in granular media — for example, oscillatory 
U EL S IU and continuous shear [1, H , horizontal Q and 
vertical vibration [H, and gravity-driven flows [To| . 
The fact that driven granular systems can achieve steady- 
states, explore configuration space, and experience fluc- 
tuations as in thermal systems, has prompted a number 
of groups to describe these flows using concepts borrowed 
from equilibrium statistical mechanics (such as effective 
temperature) [H OS 01 EI EE G1 • 

However, before a statistical mechanical treatment can 
be rigorously applied to dense granular flows, fundamen- 
tal questions about the nature of configuration space 



should be addressed. In particular, one needs to deter- 
mine how dense granular systems sample configuration 
space: Is it uniformly sampled or are some states vis- 
ited much more frequently than others? How is the sam- 
pling of configuration space affected by the strength and 
type of driving and dissipation mechanisms? In this work 
and a series of recent papers [13, El, IU, H3| , we address 
these questions with the goal of developing a comprehen- 
sive physical picture for static and slowly driven granular 
matter. 

In our previous studies, we focused on the statistical 
properties of static frictionless disk packings generated 
by slow compression without gravity [rR [l8l . fl9j or by 
gravitational deposition [20]. We have determined that 
the probability distribution for mechanically stable pack- 
ings is strongly peaked around the value typically quoted 
for the random-close packing (RCP) volume fraction, and 
explained why RCP is obtained by many compaction pro- 
tocols. We have also found that the MS-packing proba- 
bilities are highly non-uniform, contrary to the Edwards' 
equal-probability assumption [2l| that is frequently used 
in thermodynamic descriptions of granular matter [lij]. 
In addition, we have found that the probabilities become 
more non-uniform with increasing system size [r^ |. 

In the present article, we further explore the statistics 
of granular microstates and its relevance for static and 
dynamic properties of granular materials. We focus on 
slowly sheared systems at fixed zero pressure, where the 
evolution can be approximated as a sequence of transi- 
tions between MS packings. One of our novel results is 
that slowly sheared MS packings occur as continuous ge- 
ometrical families defined by the network of particle con- 
tacts. Moreover, these geometrical families are not uni- 
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formly sampled during quasistatic shear flow. We focus 
on small systems, so that we are able enumerate nearly 
all MS packings and obtain accurate packing probabili- 
ties as a function of shear strain. Since our results indi- 
cate the need for developing an alternative approach to 
the quasi-thermodynamic descriptions based on the Ed- 
wards' ideas, we advocate here a new master-equation 
framework for understanding dense granular flows. 

This paper is organized as follows. In Sec. [Hi we pro- 
vide motivation for our investigations by introducing a 
simple phenomenological master-equation model that re- 
produces qualitatively some of the key features of slowly- 
driven granular systems. In Sees. UTTI and IIV( we intro- 
duce our model system, 2D bidisperse mixtures of fric- 
tionless, purely repulsive, soft grains, and describe the 
simulation method that we employ to generate mechani- 
cally stable packings. Here, we also clearly define the set 
of distinct mechanically stable (MS) packings in terms of 
the eigenvalues of the dynamical matrix and discuss their 
symmetries. In Sec. [V] we outline our method to study 
quasistatic simple shear flow of frictionless disks in small 
2D systems at zero pressure. In Secs. lVIl we describe the 
results of the quasistatic shear simulations, with a partic- 
ular emphasis on enumerating geometrical families of MS 
packings and determining how they are sampled during 
quasistatic dynamics. In Sec. IVIll provide an outlook for 
further work on larger system sizes. The main conclu- 
sions of our studies and their relation to our long-term 
research program are discussed in Sec. I Villi In Appen- 
dices |A] |Bl and we provide details of the numerical 
simulations, the calculation of the dynamical matrix for 
the repulsive linear spring potential used in these studies, 
and the method used to distinguish 'polarizations' of MS 
packings (i.e., MS packings that differ only by reflection 
or rotation). 



II. MOTIVATION 

Granular materials are athermal — they are unable to 
thermally fluctuate and sample phase space. However, 
if a grain packing is perturbed by external forces, it can 
move through a series of configurations. The set of states 
populated by a granular system during a series of discrete 
vertical taps was characterized in Ref. 0]. As shown in 
Fig. HJa), an initially loose packing is compacted by tap- 
ping first gently, and then with greater intensity. At suf- 
ficiently large tapping intensities, it is no longer possible 
to further compact the system. However, when the tap- 
ping intensity is reduced, the packing fraction increases, 
rather than returning along the original packing fraction 
trajectory. This new curve (packing fraction vs. tap- 
ping intensity) obtained by successively decreasing and 
then increasing the tapping intensity in small steps is 
nearly reversible. A similar phenomenon has been found 
in granular media undergoing cyclic shear [f|, as shown 
in Fig.QJb). 

These experiments, which show that slowly driven 
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FIG. 1: (a) Packing density p for mm glass beads as a func- 
tion of tapping intensity P (normalized by g). The system is 
initialized in a dilute state at p = 0.59 and then subjected 
to vibrations of increasing intensity. At each P, the system 
is tapped until it achieves a steady-state p. After reaching 
r = 7, the tapping intensity is decreased until T < 1, and then 
increased again. Data is reprinted from Ref. [8|. Copyright 
1998, The American Physical Society, (b) Packing fraction cj> 
for mm glass beads as a function of shear angle 6 during cyclic 
shear. 9 is increased linearly in time from 2.7° to 10.7°, then 
decreased to 2.7°, and finally increased again to 10.7°. Data 
is reprinted from Ref. [5]. Copyright 2000, Springer- Verlag. 



granular systems appear to explore a well-defined set of 
states reversibly, have prompted a number of theoreti- 
cal studies aimed at describing compacting granular sys- 
tems using quasi-thermodynamic approaches based on 
the Edwards' ensemble {i.e. the assumption of equally 
probable microstat es) I22l.l23l.j23 . l25l |2a]. However, in 
previous work [13, fl8L Il9l. l20l|. we have shown explic- 
itly for small systems that the probability distribution 
for mechanically stable packings is highly non-uniform. 
Moreover, we have demonstrated that this feature is not 
sensitive to the packing-generation protocol and becomes 
more pronounced as the system size increases. Thus, we 
argue that the Edwards' equal-probability assumption is 
not valid and alternate theoretical approaches for slowly 
driven granular systems must be developed. 

Although several alternatives have been put forward 
, we advocate here for a master-equation ap- 
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FIG. 2: Model for quasistatic evolution of slowly driven gran- 
ular materials. Evolution of a system undergoing a long se- 
quence of quasistatic excitations (taps). The excitation am- 
plitude A is periodically ramped up and down in time t, as 
shown in the inset. The system initially evolves along the ir- 
reversible branch I, but at subsequent periods of ramping it 
moves along the reversible branch R (with slight hysteresis). 
The results were obtained by solving a set of master equations 
([3]) with model transition rates Wki (between states k and I) 
adjusted to qualitatively reproduce the packing fraction ver- 
sus tapping intensity obtained in experiments shown in Fig. [1] 
(a). 

proach for the following reasons. First, quasistatic evo- 
lution of slowly driven granular systems can be approx- 
imated as a sequence of transitions between MS pack- 
ings. Second, since the system undergoes particle re- 
arrangements as it transitions between MS packings, it 
retains little memory from one MS packing to the next, 
and successive MS packings are nearly statistically inde- 
pendent. Thus, slowly driven granular systems can be 
approximated as a Markov chain of independent transi- 
tions between MS packings and described using a master- 
equation approach. 

We have shown that the form of the microstate prob- 
ability distribution can be qualitatively reproduced by 
combining probabilities for small subsystems. Thus, we 
also advocate a novel 'bottom-up' approach in which 
large granular systems are described as collections of 
nearly independent subsystems. We view this work on 
small quasistatically sheared MS packings as laying the 
groundwork for future studies that will apply the master- 
equation approach to quantitatively describe the statis- 
tical properties of dense granular flows. 



A. Model 

To illustrate the importance of the above-mentioned 
features in capturing the irreversible and reversible 
branches in Fig. [TJ we construct a simple model in which 
a granular system is represented as a collection of N sta- 
tistically independent small subdomains. Each subdo- 
main j = 1, . . . , TV can reside in one of several microstates 
kj = 1, . . . , m. The volume of the subdomain j in state 
kj is Vi(j, kj), and these subvolumes are assumed to be 



additive 

N 

V(A) = J2Vx(j,k j ), (1) 
j=i 

where V(A) is the volume of the whole system in the state 
A = (kx, ■■ ■ , fcjy). Since the subdomains are assumed to 
be statistically independent, the joint probability distri- 
bution for the whole system is 

N 

P(A;t) = l[P j (k j ;t), (2) 

j=i 

where Pj(kj] t) is the probability distribution for the mi- 
crostates of subdomain j at time t. 

The evolution of subdomains in a system driven by an 
applied force of strength A is described by N independent 
master equations 

m 

P j (k;t i+1 )=P j (k;t i ) + J2l W kk>WPj( k ';ti)- 
fc'=i 

WvkMPjfaU)], (3) 

where Wkk' is the transition probability from state k! to 
state k. 

To qualitatively reproduce the irreversible and re- 
versible branches of states when the magnitude of the 
external forcing A is ramped up and down, we use simple 
assumptions regarding the volume distribution of indi- 
vidual subsystems and the transition probabilities. We 
assume that the volumes of individual subsystems are 
given by the expression 

Vt =r(V + Ae- K ^), (4) 

where Vo, A, and K v are constants, Xi = (i — l)/(k — 1), 
and r is a random number. Note that states with higher 
indexes possess smaller volumes. The transition rates are 
modeled by the expression 

W ij =Cb(x ij )r w e-^^ x ^ 2 , (5) 

where C is the normalization constant, b(xij) is a bias 
function 

6(x) = |6 min + 6o(l-2A) *>0, (g) 

with constants £> m i n and bo, which introduces asymmetry 
between up and down jumps, <jq determines the width 
of the Gaussian distribution, and r w is a random num- 
ber. As shown in Fig. [21 such a simple master-equation 
approach is able to qualitatively reproduce features (evo- 
lution of the system to a reversible branch of the pack- 
ing fraction) found in vertical tapping and cyclic shear 
studies of granular materials. The success of this sim- 
ple model emphasizes two key points: 1. Only minimal 
constraints on transition probabilities are required (not 
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thermodynamics) to reach steady-state and 2. An as- 
sumption of weakly interacting small subsystems may be 
able to explain macroscopic phenomena in slowly driven 
granular media. 

These results suggest a new approach to describe the 
quasistatic evolution of granular systems: a Markov pro- 
cess characterized by transfer rates between MS packing 
microstates. To further develop this approach so that it 
can yield quantitative predictions, one must determine 
(a) the types of microstates that occur in static granu- 
lar systems and (b) transition probabilities between these 
microstates when the system is slowly driven. To provide 
the necessary input for constructing quantitative descrip- 
tions, in this paper we study the structure of configura- 
tion space and transitions between microstates in small 
2D granular systems undergoing quasistatic shear strain 
at fixed zero pressure. 

III. SMALL PERIODIC GRANULAR 
PACKINGS IN SIMPLE SHEAR 

A. Bidisperse frictionless disks 

We consider 2D systems of soft, frictionless disks in- 
teracting via the pairwise additive purely repulsive linear 
spring potential 

N 

V(0= E ^(r«), (7) 

Kj = l 

Vzinj) = — (l - r lj /<7 lj ) a e (o-tj/nj - 1) , (8) 

a 

where a = 2, £ = (fx, . . . , fjv) denotes the system con- 
figuration, Ti is the position of the center of disk i, 
r ij = \?i ~ fj| is the center-to-center separation between 
disks i and j, and the sum is over distinct particle pairs. 
The strength of the spring potential ((8|) is characterized 
by the energy scale e, and the range by the average par- 
ticle diameter o~ij = (cr; + cr?) /2. The Heaviside step 
function <d(x) turns off the interaction potential when 
the particle separation is larger than a^j . 

All numerical simulations described in this paper were 
performed for 50-50 (by number) binary mixtures of large 
and small particles with diameter ratio 1.4. In such bidis- 
perse mixtures, shear-induced crystallization is inhibited 
[29( 1 : thus the system is well suited for investigations of 
quasistatic evolution of disordered granular systems. We 
focused on small systems with the number of particles in 
the range N = 4 to 20. 

To mimic the behavior of granular packings, we con- 
sider MS disk configurations at infinitesimal pressure and 
particle overlaps. As shown in our recent experimen- 
tal and numerical study, statistical properties of disks 
interacting via the repulsive linear spring potential ([8]) 
closely match properties of plastic and steel disks in a 
system where frictional forces have been relaxed using 
high-frequency, low-amplitude vibrations (20| . 
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FIG. 3: Schematic of shear-periodic boundary conditions, (a) 
MS packing with N = 6 particles confined to a L x L box with 
shear strain 7 = Ax/L — 0. (b) MS packing in a unit cell 
with 7 = 0.2. Note that at arbitrary 7, a given particle in the 
primary cell is not directly above (or below) its image, (c) 
MS packing at 7 = 0.8, which shows that configurations at 7 
and 1 — 7 are related by an inversion about the vertical axis, 
(d) MS packing in (a) at 7 = 1. At unit shear strain, shear- 
periodic boundary conditions are identical to standard peri- 
odic boundary conditions [3^]. Thus shear-periodic boundary 
conditions have unit period. In all panels, the shaded parti- 
cles indicate a given particle in the primary cell and its image 
in a neighboring cell. 

B. Shear-periodic boundary conditions 

In our simulation studies, the particles are confined to 
a L x L periodic box, as illustrated in Fig. [3] The unit 
cell can either be a square [cf. Fig. 03a)] , or it can be de- 
formed in the x direction [cf., Figs.^b) and (c)]. These 
shear-periodic boundary conditions allow us to generate 
an ensemble of anisotropic granular packings as a func- 
tion of the shear strain 7 = Ax/L, where Ax is the 
horizontal shift of the top image cells relative to the bot- 
tom image cells. Moreover, simulations of systems with 
gradually changing strain [3(| Hl| enable us to study qua- 
sistatic evolution of a granular packing under shear. 

Note that shear-periodic boundary conditions are iden- 
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FIG. 4: Schematic of the compression/decompression pro- 
tocol to create MS packings. In panels (a)-(c), we show a 
schematic of the energy landscape V(£) in the vicinity of the 
static granular packings (at point £o in configuration space) 
in panels (d)-(f). If the system exists in a non-overlapped 
configuration (panel (f)) with gaps between particles and 
zero energy per particle (V = 0), it will be compressed fol- 
lowed by energy minimization. If the system exists in an 
overlapped configuration (panel (d)) at a local energy min- 
imum with V > 0, it will be decompressed followed by en- 
ergy minimization. When the system switches between the 
cases displayed in panels (d) and (f), the size of the com- 
pression/decompression increment is decreased. The process 
stops when the system exists in a static granular packing at 
a local energy minimum that is infinitesimally above zero. 
This schematic is shown for shear strain 7 = 0, but a similar 
process occurs for each 7. 

tical at 7 and 1 — 7 as shown in Fig. [3l Also, when 
reflection symmetry is taken into account, it is clear 
that we only need to consider shear strains in the range 
7 = [0, 0.5] to generate static packings at arbitrary shear 
strains. However, in the case of continuous quasistatic 
shear flow, evolution of the system over multiple shear 
strain units must be investigated to capture the full dy- 
namics. 











FIG. 5: (I) A typical MS packing for N = 6 particles at inte- 
ger shear strain. The configurations in panels (2)-(8) are ob- 
tained by applying the seven possible reflections and rotations 
in two spatial dimensions consistent with periodic boundary 
conditions in an undeformed square cell. (See Appendix [C]) 
Configurations in opposing columns are related by a rotation 
R.3 by n about an axis coming out of the page, (bottom) 
This schematic shows that configurations undergoing simple 
planar shear in 2D are invariant with respect to R3. 



TABLE I: Number of distinct MS packings at integer shear 
strains when we treat all polarizations as the same N a and 
when we distinguish among different polarizations Nf as a 
function of system size N. The third column gives the ratio 
Nf/N B . Data for N s is obtained from Ref. \v%. 



_/V Ns Nf Nf/Ns 

4 3 6 2.00 

6 20 68 3.40 

8 165 612 3.71 

10 1618 6378 3.94 

12 23460 91860 3.92 



IV. ENUMERATION OF MS PACKINGS AT 
ARBITRARY SHEAR STRAIN 

A. Packing-generation protocol 

Zero-pressure MS packings at a given shear strain are 
generated using the compression/decompression packing- 
generation protocol used in our previous studies of un- 
sheared MS packings We briefly outline the pro- 

cedure below for completeness. We begin the packing- 
generation process by choosing random initial particle 
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FIG. 6: The number of distinct mechanically stable packings 
at integer shear strain when we treat all polarizations the 
same N 3 (circles) and when we account for different polar- 
izations Nf (squares). The solid and long-dashed lines have 
slope « 1.2. 

positions within the simulation cell at packing fraction 
d>a = 0.50 (which is well below the minimum packing 
fraction at which frictionless MS packings occur in 2D). 
We then successively increase or decrease the diameters 
of the grains, with each compression or decompression 
step followed by conjugate gradient minimization [32j of 
the total energy in ([7]). 

As illustrated in Fig. [4] the system is decompressed 
when the total energy (|7J) at a local minimum is 
nonzero — i.e., there are finite particle overlaps [cf., 
Figs. H] (a) and (d)]. If the potential energy of the sys- 
tem is zero [Fig. Q|c)] and gaps exist between particles 
[Fig. Hff)], the system is compressed. The increment by 
which the packing fraction <fi is changed at each com- 
pression or decompression step is gradually decreased. 
Numerical details of the algorithm are provided in Ap- 
pendix [XJ 

In the final state of the packing-generation process, the 
potential energy vanishes [Fig. 1Kb)] , but any change of 
the relative particle positions (excluding rattler particles, 
which can be moved without causing particle overlap) re- 
sults in an increase in the potential energy. Thus, the 
final state is a mechanically stable configuration (or col- 
lectively jammed state [33| ) at jamming packing fraction 
At each fixed 7, MS packings form a discrete set in 
configuration space. The packing-generation process is 
repeated more than 10 6 times for at least 100 uniformly 
spaced shear strain values in the interval < 7 < 0.5. A 
large number of independent trials is required to enumer- 
ate nearly all MS packings because the MS probability 
distribution varies by many orders of magnitude. 




FIG. 7: The number of distinct mechanically stable packings 
Nf (7) as a function of shear strain 7 for N = 4 (long-dashed) , 
6 (solid), and 10 (dotted). Nf^y) is normalized by the average 
number of MS packings over all shear strain (Nf^)}. The 
sampling interval is A7 = 10~ 2 . As shown in Fig. [3] the 
same set of MS packings occur at 7 and 1 — 7. 

it is very rare that two distinct MS packings have the 
same 4>j. Thus, it is often convenient to characterize MS 
packings by cf>j. However, in our detailed investigations 
of the quasistatic evolution of sheared granular packings, 
a more precise classification of MS packings is necessary. 

To determine the set of distinct MS packings, we use 
the eigenvalues and eigenvectors of the dynamical matrix 

SI 



dV 



(9) 



£0 



where £ m is the mth component of the configuration vec- 
tor £ and £0 gives the configuration of the reference MS 
packing. Since rattler particles do not contribute to the 
stability of the system, we consider only the dynamical 
matrix for the mechanical backbone of the packing. This 
matrix has dAf rows and columns, where d = 2 is the 
spatial dimension, Af = N — N r is the number of parti- 
cles in the mechanical backbone, and N r is the number 
of rattler particles. 

Since the dynamical matrix is symmetric, it has dM 
real eigenvalues {rrii}, d of which are zero due to trans- 
lational invariance of the system. In a mechanically sta- 
ble disk packing, no collective particle displacements are 
possible without creating an overlapping configuration; 
therefore the dynamical matrix for MS packings has ex- 
actly d(Af — 1) positive eigenvalues [35(]. We limit the 
results below to mechanically stable packings. 



B. Classification of granular packings 

1. Dynamical matrix 

With our precise measurement of the jamming pack- 
ing fraction <f>j to within 10~ 8 of the jamming point, 



2. Polarizations 

According to our previous investigations [13, HH , dis- 
tinct MS packings always have distinct sets of eigen- 
values {mi}, except when packings can be related to 
each other by reflection or rotation (35|. When we treat 
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(b) 




FIG. 8: Schematic of the evolution of the energy landscape 
during quasistatic shear at fixed zero pressure from shear 
strain 7 to y+5j. In (a), the system evolves continuously from 
the local minimum at shear strain 7 (1) to the one at 7 + £7 
(4) because there are no particle rearrangement events dur- 
ing the shear strain interval. In contrast, in (b) we show that 
when the system undergoes a particle rearrangement event, 
it will reside in a different energy minimum at 7 + £7 com- 
pared to the one at 7. Snapshots of the static packings at 
shear strain 7 (gray) and 7 + 57 (black) are superimposed in 
(c) and (d), which correspond to the potential energy land- 
scape dynamics in (a) and (b), respectively. In (d), three of 
the original contacts are removed and four new contacts are 
generated as a result of the particle rearrangements. 




y 

FIG. 9: The jammed packing fraction <pj versus shear strain 
7 during quasistatic shear at zero pressure for N = 10. The 
dashed line shows the initial transient and the solid line high- 
lights the periodic behavior with unit period that begins near 
7t ~ 2.5. The inset shows <j>j for the same system in the main 
figure except only at integer shear strains. 



Distinct MS packings at integer and non-integer 
strains 



different 'polarizations' associated with each symmetry 
transformation as the same state, distinct MS packings 
can be classified according to the lists of eigenvalues of 
their dynamical matrices. This approach was adopted in 
Refs. [T3, Ell], where we considered only static, isotropic 
particle configurations at 7 = 0. The eight equivalent 
polarizations for a MS packing at 7 = are shown in 
Fig. H 

In contrast, in our present study we consider continu- 
ous shear deformations of the system. After an isotropic 
unit cell is deformed, different polarizations of a given 
state can be transformed into distinct MS packings (i.e., 
packings distinguishable by the lists of eigenvalues {m.;}). 
For example, all polarizations shown in the left (right) 
column of Fig. O after a step strain are transformed into 
non-equivalent configurations. Our present classification 
scheme for MS packings takes this effect into account. 

Accordingly, we treat states that differ by reflection 
or rotation by the angle a = it/ 2 or 3ir/2 as distinct 
MS packings. However, due to symmetry of shear flow 
(cf. bottom of Fig. [5]), the states that are related by a 
rotation by the angle a = it deform in equivalent way. 
Therefore, such states are treated as equivalent states. 
Further details regarding polarizations of MS packings, 
including the procedure we use to distinguish polariza- 
tions, are discussed in Appendix [Cl 



MS packings at integer shear strains have equivalent 
boundary conditions to those in standard square periodic 
cells (as shown in Fig. [3]). Thus, we can use our previous 
extensive calculations of MS packings generated in square 
cells with periodic boundary conditions [l7l ]. 



In TableQ]we show the number of distinct MS packings 
at integer shear strains N s when we treat all polariza- 
tions as the same (data from Ref. [13]) and Nf when we 
account for different polarizations (as described above). 
The ratio Nf /N s is smaller than 4 because of reflection or 
rotation symmetry of some MS packings. As depicted in 
Fig. [HI both N s and Nf grow exponentially with system 
size. 



The number of distinct MS packings Nf(-f) (includ- 
ing polarizations) versus shear strain is shown in Fig. 
for several system sizes. The results show that (1) the 
maximum number of distinct packings occurs at integer 
and half-integer shear strains, and (2) there is a notice- 
able dip in Nf at 7* w 0.2. However, as the system size 
increases, the number of distinct MS packings becomes 
uniform as a function of shear strain. 
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TABLE II: The maximum transient shear strain 7t nax and pe- 
riod T of each distinct periodic orbit in integer strain units 
for several system sizes N. (7t" ax is rounded to the nearest 
integer shear strain.) In the fourth column, we list the av- 
erage jammed packing fraction (</>./) of the MS packings that 
are dynamically accessible at large shear strain for each each 



periodic orbit. 


For systems with multiple 


periodic orbits, we 


list 7t max , T, and (cf> 
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TABLE III: Catalog 
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V. QUASISTATIC SHEAR FLOW AT ZERO 
PRESSURE 

A. System dynamics 

1. Quasistatic shear-strain steps 

To mimic quasistatic evolution of a frictionless gran- 
ular packing, we apply a sequence of successive shear- 
strain steps of size <5 7 <C 1 to a system of bidisperse disks 
with shear-periodic boundary conditions. Each step of 
the protocol consists of (1) shifting the ^-coordinate of 
the particles, 

Xi^Xi+Sjyt, (10) 

in conjunction with the corresponding deformation of 
the unit cell; and (2) the compression/decompression 
packing-generation process (described in Sec. IIV Aj) to 
achieve a zero-pressure MS configuration with infinitesi- 
mal particle overlaps. 

This procedure generates quasistatic shear flow at con- 
stant zero pressure (but varying packing fraction) , which 
is an appropriate description of slowly sheared granu- 
lar matter, where particle overlaps are always minimal. 
During the evolution, the system constantly expands and 
contracts to remain at the onset of jamming with pack- 
ing fraction <j)j that depends on 7 . We note that our 
procedure is distinct from the quasistatic evolution used 
in recent investigations of sheared glass-forming liquids 
[HI, [37j . where a constant-volume ensemble was imple- 
mented. 



2. Particle rearrangements 

During a single shear-strain step 5j, particle positions 
can either change continuously or a sudden particle re- 
arrangement can occur. A continuous evolution step is 
schematically depicted in the left panels of Fig. [51 and the 
right panels illustrate a strain step that yields a particle 
rearrangement . 

In both cases, the system initially resides in the lo- 
cal minimum (point 1) in the energy landscape. Since 
this minimum corresponds to a MS packing, its energy 
is infinitesimal. When the affine transformation (|10p is 
applied, the system moves to point 2 in the energy land- 
scape. In the case of continuous evolution, depicted in 
Fig. [SJa), the energy minimization (point 3) and subse- 
quent decompression (point 4) move the system back to 
the neighborhood of the initial packing. As illustrated 
in [8] (c), at the initial and final points (1) and (4), the 
topology of contact networks is unchanged. 

In contrast, when the energy minimization that fol- 
lows an affine shear-strain step drives the system into a 
zero-energy region corresponding to an unjammed pack- 
ing (i.e. point 3 in Fig. [SJb)), upon subsequent com- 
paction the system is driven to a MS packing (point 4) 
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FIG. 10: This 'tree' diagram for JV = 6 shows the evolution 
at integer shear strains of the system undergoing quasistatic 
shear flow at zero pressure. The system is initialized at all 68 
possible MS packings using the same indexes i as in Table HTT1 
The shear strain axis has been shifted so that the periodic 
orbits for all initial MS packings begin at 7 = 2. There is 
only one MS packing (67) that is dynamically accessible at 
large shear strain. 

with a different contact network than the initial packing 
(point 1), as illustrated in Fig.[5[d). 

During the continuous portion of the quasistatic shear 
strain evolution, the MS packings that are visited do 
not depend on the energy minimization method (e.g., 
energy relaxation via dissipative forces versus the con- 
jugate gradient algorithm [171 ]) or parameters related to 
the compaction and decompression processes. Dynami- 
cal features do not influence the MS packings that are 
obtained along the continuous region because the system 
remains in the basin of the original mechanically stable 
packing. However, when a shear strain step leads to par- 
ticle rearrangements, the system is taken to a new region 
of the energy landscape, and the energy minimization 
scheme, compaction and decompression rates, and even 
the location of rattler particles can influence the final MS 
packing. In small systems (N < 16), we find that the 
dynamics is weakly sensitive to these features, whereas 
in larger systems noise and protocol dependence play an 
important role in determining steady-state MS packing 
probabilities. In Sec. IVII1 we will describe future work 
in which we will tune the dynamics to determine its in- 
fluence on the transition rates among MS packings. 



B. Deterministic and contracting evolution of 
small frictionless MS packings 

Our key results regarding quasistatic evolution of small 
granular systems are summarized in Figs . [91 and 1101 and in 
Tables [TTI and IIII1 The results show that (1) the evolution 
of small systems is deterministic; (2) after transient evo- 
lution, the system becomes locked into a periodic orbit; 



and (3) the evolution in configuration space is strongly 
contracting, in the sense that each unit strain leads to a 
significant reduction of the number of dynamically acces- 
sible MS packings. 

Deterministic evolution The deterministic character 
of the system evolution over continuous portions of the 
trajectories can be justified using arguments based on 
the continuity of the energy landscape, similar to those 
illustrated in Fig. [5] (a). The exceptions are bifurcation 
points, which will be described in Sec. IVIIl Indeter- 
minacy can also occur due to the presence of rattlers 
or random particle motions in unjammed configurations. 
For systems with N < 16 evolved according to the al- 
gorithm described in Sec. IV A 11 the observed evolution 
was always completely deterministic. Random evolution 
in larger systems and systems with random noise are dis- 
cussed in Sec. IVIIl 

Periodic orbits In Fig. we track the evolution of 
the jammed packing fraction (f>j as a function of shear 
strain 7 during zero-pressure quasistatic shear flow after 
the system is initialized in one of the MS packings at 
7 = for N — 10 particles. Fig. [5] shows the complete 
trajectory and the inset shows <pj only at integer strains. 
We observe that after a short transient of approximately 
two units of strain, the system becomes locked into a pe- 
riodic orbit (with unit period T = 1). Similar behavior is 
observed for systems with N = 4 to 14 particles, as sum- 
marized in Table [TTI Although we have a limited range 
of system sizes for which a complete analysis of the qua- 
sistatic evolution has been preformed, the results show 
that both the transient time jt and period T increase 
somewhat with system size. (Note that several periodic 
orbits for N = 12 and 14 particles have anomalously large 
transients, cf. Sec. IVIIO 

Contracting evolution The contracting character of 
the quasistatic shear flow is illustrated in the 'tree' di- 
agram in Fig. 1101 This digram represents the complete 
set of trajectories, shown for integer strain, for a sys- 
tem oi N — 6 particles. For this system size there are 
Nf — 68 possible MS packings at any given integer shear 
strain, yet only one of them is dynamically accessible in 
the large shear strain limit. 

The contraction of the set of dynamically accessible 
states as the system evolves results from the determinis- 
tic and irreversible character of quasistatic evolution. As 
illustrated in Fig. [TJJl from each MS packing at an integer 
strain 7, the system transitions to a unique MS packing 
at the strain 7 + 1; however, more than one state can 
transition to a given state. A complete catalog of tran- 
sitions at integer shear strains for the system of A = 6 
particles is given in Table IIIII From this table, a com- 
plete list of trajectories of the whole evolution process 
can be constructed without further simulations. 

As shown in Table [TTJ the number of periodic orbits 
to which the system contracts increases with the system 
size. However, for all systems studied, the number of 
orbits and the number of states visited at long times is 
small compared to the total number of MS packings. 
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FIG. 11: Schematic of the evolution of the jamming packing 
fraction <j>j during the shear strain interval 7 S to 73+1. In 
(a) the particle contact network does not change from "f s to 
7s+i, while there is a discontinuity in cj>j in (b) and in the 
derivative of (j>j in (c) at 7*, which signal a change in the 
contact network. 

The results in the last column of Table [H] give the 
average jammed packing fraction (<pj) of the MS packings 
sampled for a given periodic orbit. The MS packings that 
occur in periodic orbits typically have the highest packing 
fractions out of the entire distribution. 



VI. FAMILY STRUCTURE OF THE SET OF MS 
PACKINGS 

An analysis of the deterministic evolution of the sys- 
tem monitored at integer values of strain is sufficient to 
conclude that the dynamics is contracting and periodic 
at large shear strains. However, it might be puzzling that 
both the transient shear strain and period are so short, 
i.e., much shorter than the number of MS packings (cf., 
Tables HI and Hi]) . To shed light on this behavior we will 
now analyze how the sets of MS packings at different 
values of shear strain are connected. 



TABLE IV: Statistics for complete geometrical family maps. 
Total number of distinct geometrical families Nf, average 
number of distinct MS packings {Nf (7)), average packing 
fraction {4>j(j)), average family length f, and average fam- 
ily second derivative C for several system sizes N 38|. 





N 


N f 


(JV?( 7 )> 


(Ml)) 




C 


Y 


4 


15 


4 


0.788 


0.18 


1.7 




6 


334 


47 


0.739 


0.06 


5.6 




10 


34822 


2896 


0.753 


0.03 


4.9 



family to another) is accomplished either by jumps (dis- 
continuities in <j)j) or kinks (discontinuities in the deriva- 
tive of 4>j with respect to 7) as shown in Figs. QT] (b) 
and (c), respectively. As will be discussed in Sec. IVIB] 
a discontinuity in <pj corresponds to a system instability 
and a change in the contact network with finite particle 
displacements; whereas a discontinuity in the derivative 
of 4>j corresponds to a change of the contact network 
without finite particle displacements. 

To construct the complete map of distinct geometrical 
families for all shear strain, we can simply link the equiv- 
alent geometrical families at the shear strain endpoints 
7 S , or terminate the family if it has no counterpart. Be- 
cause of the contracting dynamics described in Sec. IV B| 
it is important to use sufficiently small shear strain in- 
tervals so that we have enough shear strain resolution to 
capture small families. (See Appendix [5] for additional 
details for constructing geometrical families.) 



Complete map of MS packings for continuous 
shear strain 



A. Construction of continuous geometrical families 
of MS packings 

As described in Sec. IIVB1 we can identify the distinct 
MS packings at a given shear strain using the eigenvalues 
of the dynamical matrix. However, this method will not 
work for comparing MS packings at different shear strains 
since the eigenvalues vary continuously with 7. 

To study the relationship between distinct MS pack- 
ings at different shear strains, we divide the shear strain 
region into small intervals 7 s +i — 7s- For each distinct 
MS packing at 7 S , we monitor the particle contact net- 
work as the system evolves toward 7 s +i (and 7 s -i). In 
Fig. QT] (a), we show the continuous evolution of <j)j be- 
tween 7 S and 7 s +i, which implies that there are no rear- 
rangement events and no changes in the particle contact 
network during this interval. Thus, the continuous evo- 
lution of 4>j identifies a portion of a geometrical family of 
MS packings all with the same particle contact network 
that exist over a continuous shear strain interval. 

In our simulations, we find that changes in the network 
of particle contacts (i.e. switches from one geometrical 



We find that a particularly simple, pictorial method 
to distinguish geometrical families is by monitoring the 
jamming packing fraction cf>j as a function of shear strain 
(instead of a complete representation of the particle con- 
tact network). The complete map of MS packings — <pj 
for all distinct MS packings at all shear strains — is dis- 
played in Fig. [T^] for N — 6. The structure of the map is 
quite complex even for TV = 6, and it possesses a number 
of distinctive features. 

First, the map is composed of a finite number (Nf — 
334 for N = 6) of curved concave-up segments each of 
which corresponds to a distinct geometrical family of MS 
packings. Second, the parabola-like curves either end dis- 
continuously [cf., point 2 in the blowup shown in Fig. 1121 
(b)] or form a kink [points 4 and 5 in Fig. [121(b)]. Third, 
the curves significantly vary in length, and they are not 
symmetric about the apexes (which are distributed over 
the full range of 7 [3!| ) . 

We find that a given parabola-like curve ends when 
either the contact network of that particular continuous 
region becomes unstable and the system jumps to a new 
one or the contact network merges with another network 
to form a kink. Examples of contact networks corre- 
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FIG. 12: (Color online) (a) Complete geometrical family 
map — jamming packing fraction (j>j("i) as a function of shear 
strain 7 — for N = 6 (blue lines) . (b) Close-up of boxed region 
in (a). The orange (black) lines indicate increasing (decreas- 
ing) shear strain evolution, (c) Particle configurations at five 
points during evolution within the family map are also shown. 
Solid lines connecting particle centers represent particle con- 
tacts; each distinct network is given a different grayscale. 
Contacts denoted by thick lines in panels 4 and 5 are either 
gained or lost as the system evolves from configuration 4 to 
5. 

sponding to characteristic points on the family map are 
shown in Fig. [T2] (c). 

Figure [TUl shows that the family-length distribution is 
exponential (with the characteristic strain of £ ~ 0.06 for 
N=6 |4(|) and. the distribution of second derivatives pos- 
sesses a strong peak. We speculate that the decay length 
£ is related to the average yield strain in frictionless MS 
packings. Thus, it is important to study £ as a function 
of system size and this direction will be pursued in future 
work. 

The jammed packing fraction <pj has parabolic-like de- 
pendence on 7 because we consider jammed disk pack- 
ings. The general feature of these packings is that 
as shear strain first increases the system must dilate 
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FIG. 13: (a) Probability distribution P(Lj) for the length (in 
units of strain) of geometrical families for N = 6. The solid 
line represents exponential decay with decay length £ ~ 0.06. 
(b) Distribution P(C) of second derivatives of 0,7(7) (with 
respect to shear strain 7) for the distinct geometrical families 
for N = 6. The family average is C ~ 5.6. 

(packing fraction decreases) as particles climb over each 
other [4l|] . However, beyond a critical shear strain the 
system must compact to maintain particle contacts. 

Since the shape of the family map 4>j — 4>j("f) is en- 
tirely determined by the geometry of the particle contacts 
at zero pressure, this shape is independent the detailed 
form of the finite-range potential l[8|). In particular, it 
does not depend on the power a defining the interparti- 
cle elastic forces. 

We note that recent simulations have found that pack- 
ings of ellipsoidal particles possess a large number of 
low-energy vibrational modes, which are quartically, not 
quadratically, stabilized [52, SH, |44| . Thus, we predict 
non-parabolic dependence of </>j(7) for small quasistati- 
cally sheared packings of ellipsoidal particles. 

In Table IIVI we summarize our results by compiling 
statistics of the complete geometrical family maps for 
several system sizes N. We provide the total number 
of distinct geometrical families Nf, average number of 
distinct MS packings (iVf (7)), average packing fraction 
(</>j(7)), average family length £, and average family sec- 
ond derivative C (with respect to 7), where (.) indicates 
an average over shear strain 7. 



C. Quasistatic evolution 

Since we have now mapped all of the MS packings over 
the full range of shear strain, we can now relate key fea- 
tures of the quasistatic evolution at zero pressure to the 
family structure of the MS packings. We recall from Sec. 
[V]that the evolution is strongly contracting, non-ergodic, 
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FIG. 14: The complete geometrical family map as shown in 
Fig. [12] as a function of shear strain 7 for TV = 6 (gray lines). 
The black line shows the evolution of <j>j during quasistatic 
shear strain starting from a single MS packing with cj>j = 
0.742 at 7 = 0. 
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FIG. 15: The complete geometrical family map as shown in 
Fig. [12] as a function of shear strain 7 for N = 6 (gray lines). 
In (a) and (b), we show the evolution of 4>j{.l) under qua- 
sistatic shear strain starting from all distinct MS packings at 
7 = and 0.23, respectively (black lines). The number of 
distinct MS packings at 7 = 0.23 is smaller than the number 
at by a factor of « 3. 
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FIG. 16: Magnifications of the complete family maps and 
quasistatic shear evolution at 7 — 70 for initial shear strains 
(a) 70 = and (b) 0.23 from Figs. fa) and (b). 

and after a short transient, the system settles on a peri- 
odic trajectory. 



1. Topological changes and transitions between geometrical 
families of MS packings 

We first examine a detailed snapshot of the evolution 
depicted by orange line in Fig. [12] (b). Characteristic 
points along the trajectory are indicated by the circles 
marked 1-5, and the corresponding contact networks of 
the evolving MS packing are shown in (c). 

The evolution starts at point 1, which represents one 
of the MS packings in a system with an unstrained unit 
cell. When the strain is gradually increased, the MS 
packing evolves continuously from point 1 to 2 along the 
geometrical family that includes the initial point. As 
shown in (c), there is no topological change of the contact 
network during the continuous part of the trajectory. 

When the strain is increased beyond point 2 (where 
the family ends), the system discontinuously transitions 
to another branch of MS packings. During the transi- 
tion, the packing fraction <f>j discontinuously increases, 
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FIG. 17: The number N v of distinct MS packings visited 
during quasistatic shear (normalized by Ng(yo)) versus 7 — 70 
for N = 6 obtained obtained by summing over all of the 
possible MS packings initialized at 70 for 70 = (dashed 
line) and 0.23 (solid line). 

particles undergo finite displacements, and a new contact 
network is formed. Due to the condition that the system 
must be mechanically stable at zero pressure, only jumps 
up in the packing fraction are allowed. 

With further increases in shear strain, the evolution 
remains continuous until point 4, where the system en- 
counters a kink in (jud) (i-e., discontinuity in the deriva- 
tive of (f>j(j)). The kink signals a change in the particle 
contact network, with no finite particle displacements. 
Such a change occurs when a new interparticle contact 
is formed and another is broken to prevent the system 
from being overconstrained, as shown in panels 4 and 
5 in Fig. [12] (c). A change in the direction of particle 
motion also accompanies a kink in <j)j. Note that kinks 
provide an important mechanism by which the <j)j of a 
MS packing can decrease during quasistatic shear at zero 
pressure. Continued increases of the shear strain beyond 
point 5 give rise to continuous evolution of 4>j{l) until 
the next discontinuity. Only two types of discontinuities 
in cj)j — jumps and kinks — have been observed. 

We note that jumps of the system to new contact net- 
works at terminal points of continuous families provide 
a mechanism for hysteresis and irreversibility as found 
in recent experiments of cyclically sheared granular sys- 
tems [5]. If the applied shear strain is reversed at point 
3 in Fig. [T2] (b), the system will evolve along the black- 
highlighted region of 4 > j('j)j n °t the original one sampled 
during the increasing applied shear strain. 

For small systems (N < 16) undergoing quasistatic 
shear strain at zero pressure (using conjugate gradient 
energy minimization), we have found that the process of 
jumping from an old to a new MS packing family is de- 
terministic, but further study is required to determine to 
what extent it depends on the packing generation proto- 
col. Below, we will show that bifurcations in configura- 
tion space can affect the destination of MS packings dur- 
ing quasistatic shear strain; this topic will be discussed 
m more detail in Sec. EH 



2. Periodic orbits and contracting evolution 

A representative trajectory </>j — <frj(j) over several 
strain units is shown in Fig. 1141 The system is initialized 
in one of the MS packings (<fij = 0.742) at 7 = 0, and we 
plot <fij as a function 7 overlayed on the complete family 
map for N — 6. The trajectory exhibits continuous parts 
separated by kinks and jumps, similar to those described 
in Sec. IVI C 11 After a short transient evolution for 7 < 
7 t w 1.25, the trajectory becomes periodic with period 
1, consistent with the results discussed in Sec. IV Bl 

The behavior shown in Fig.Q3]is typical for small pack- 
ings under quasistatic shear strain. In Fig. [TS] (a), we 
show the evolution of 4>j for all Nf(j = 0) = 68 MS 
packings at initial shear strain 70 = 0. After a short 
initial transient 7 « 0.4, most of the initial conditions 
have converged onto one of several persistent trajecto- 
ries. After 74 w 2.25, all trajectories have converged onto 
a periodic orbit with unit period. 

Figure Qj)] (b) shows the corresponding results for all 
NP(jo = 0.23) = 26 MS packings beginning at 70 ~ 
0.23 (which is a region of shear strain at which there 
are the smallest number of MS packings according to the 
results depicted in Fig. [7]). Qualitatively, the picture is 
similar to that in (a): there is rapid contraction, several 
persistent trajectories, and then collapse onto a single 
periodic orbit. The two ensembles with 70 = and 0.23 
sample different MS packings during the transient, but 
evolve to the same periodic orbit found in (a). Thus, 
these systems sample only a small fraction of the possible 
geometrical families in the large shear strain limit. 

To determine the source of the rapid contraction of the 
number of dynamically accessible states, we magnify the 
region of small shear strains 7 — 70 in Fig. [TfJ] The results 
indicate that the contraction occurs when more than one 
trajectory jumps to the same branch. The initial con- 
traction happens quickly as depicted in Fig. II 71 although 
there is some dependence of the contraction on 70. 

An examination of the closeup in Fig. [TH] suggests that 
the source of the quick contraction and transient to the 
final periodic orbit is threefold: (i) there is a large num- 
ber of short families (cf., Fig. [TB"]) ; (ii) jumps always 
occur up in packing fraction, which makes low-packing- 
fraction MS packings dynamically inaccessible; and (iii) 
some families collect significantly more MS packings per 
unit length than the others, as discussed below. 

To estimate the 'basins of attraction' for each geo- 
metrical family, we calculated the number of jumps that 
each geometrical family collects during quasistatic shear 
strain. In Fig. 1181 we show the number of trajectories 
Nh that jump to a particular geometrical family during 
quasistatic shear as a function of the average jammed 
packing fraction (f>j of the geometrical family for N = 6. 
Nh is normalized by the family length to account for the 
fact that longer families can collect more trajectories. As 
shown in Fig. [15] we find a general trend that families at 
higher packing fractions collect relatively more of the tra- 
jectories. 
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FIG. 18: The number of trajectories Nh that jumped to a 
particular geometrical family (normalized by the geometri- 
cal family length) during quasistatic shear versus the average 
packing fraction <j>j of the given geometrical family. Jumps 
were collected over narrow strain intervals A7 = 10~ 2 for 
TV = 6. The solid line indicates an average over bins in pack- 
ing fraction with size 0.02. 

VII. TRAJECTORY SPLITTING: 
BREAKDOWN OF DETERMINISM 




FIG. 20: A schematic of the energy landscape V(£), which 
has a flat region near local energy minima (MS packings) a 
and b. For the quasistatic shear flow simulations, numeri- 
cal precision and specific details of the energy minimization 
scheme will influence whether the system reaches point a' or 
b' in the landscape, which in turn influences the likelihood of 
the system residing in MS packing a or b. 



In the previous Sees. IVl and I VII we described our re- 
sults for quasistatic shear flow at zero pressure for small 
systems TV < 14, which displayed deterministic and con- 
tracting dynamics. In particular, we found that when a 
small MS packing undergoes a jump discontinuity at a 
given 7, it makes a transition to a uniquely determined 
MS packing. Thus, these systems lock into periodic or- 
bits with typically small periods and the number of MS 
packings that are sampled at large shear strain is a small 
fraction of all possible MS packings. 
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FIG. 19: The jammed packing fraction (j>j versus shear strain 
7 during quasistatic shear flow at zero pressure beginning 
from an initially unstrained MS packing for TV = 20. Af- 
ter a short initial transient (74 ~ 1), the system locks into 
an apparent periodic orbit with large period T > 100. The 
first cycle is denoted using long-dashed and solid black lines. 
However, during the initial phase of the second cycle (long- 
dashed line), the trajectory begins to deviate from that in the 
first cycle near 70 = 137.6 (vertical dashed line). The non- 
repeating part of the second cycle is depicted using a solid 
gray line. In the inset, we overlay cycles 1 (dotted) and 2 
(solid) near 70. 



In contrast to the results found for smaller systems, 
we find that for slightly larger systems (TV > 16) there 
is a dramatic increase in the period. In addition, the 
deterministic behavior appears to break down, i.e. we 
can not predict with unit probability the transitions from 
one MS packing to another. 

The breakdown of determinism is shown for TV = 20 
in Fig. 1101 where we plot the evolution of <f>j during qua- 
sistatic shear at zero pressure. After a short transient 
strain (74 ~ 1), the system falls onto an apparent pe- 
riodic orbit with large period T > 100. However, as 
shown in the inset to Fig. [T9l during the second cycle at 
70 ~ 137.6, the trajectory of the first and second cycles 
begin to deviate. 

There are several possible mechanisms for the intro- 
duction of stochasticity and sensitivity to initial condi- 
tions in these systems, which include bifurcations in con- 
figuration space caused by local symmetries of the MS 
packings and noise (or numerical error) from the packing- 
generation protocol [45( | . 

An example of a bifurcation in the energy landscape 
is shown in Fig. [501 This region of the landscape is ex- 
tremely flat, and thus the state of the system will depend 
on the numerical precision and specific details of the en- 
ergy minimization scheme. For example, if the energy 
minimization stops at point a' in Fig. 1201 it will likely 
move toward MS packing a during the packing-generation 
process, while if it stops at point b', the system may pro- 
ceed to MS packing b. 

The configurational view of the trajectory splitting 
mechanism in Fig. [IT)] is demonstrated in Fig. [3TJ where 
we show the configuration before and after the splitting 
event at 7 = 70 . There are subtle differences in the posi- 
tion and interparticle contacts of a single particle in the 
central region of the cell at strains 70 — T (gray outline) 
and 70 (black outline). This small change in the contact 
network, which occurs in a flat region of the energy land- 
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FIG. 21: Particle configurations immediately (a) before and 

(b) after the trajectory splitting event in Fig. 1191 In (b) and 

(c) , the configurations at shear strains 70 and 70 — T are over- 
layed. The only difference between the two configurations in 
(b) is the position and resulting interparticle contacts between 
the central and neighboring particles in the dashed box. The 
particle outline and bonds are black (gray) for the configura- 
tion at 70 (70 — T). Panel (c) is a magnified version of the 
dashed box in (b). 

scape, leads to large differences at subsequent values of 
shear strain. In this case, the cause of the flatness in the 
energy landscape stems from the fact that the directions 
fij for different contacting particles j of a central particle 
i are nearly collinear. 

Similarly, if we add random displacements during 
quasistatic shear, we can create bifurcating trajectories 
0,7(7). In Fig. [HI we show the evolution of ^7(7) for sys- 
tems with and without added Gaussian random displace- 
ments with both systems initialized with the same MS 
packing at 7 = 0. We find that the trajectory with noise 
deviates from the original trajectory at 7 « 0.2. Thus, 
noise is able to compete with the contracting mechanism 
to increase the fraction of dynamically sampled MS pack- 
ings. In future studies, we will determine the fraction of 
MS packings visited as a function of the noise amplitude 
and system size. 

Rattler particles can also lead to sensitivity to initial 
conditions. Small changes in the location of the rattler 
particle even in MS packings that otherwise have the 
same network of particle contacts can lead to large differ- 
ences in subsequent contact networks if the rattler joins 
the connected network in different locations. This effect 
occurs because there is no energetic incentive for rattler 
particles to be located in any particular location within 
the confining void region as long as it does not overlap 
another particle. The influence of rattler particles on 
transition rates will also be assessed in future studies by 
varying the noise amplitude. 



VIII. CONCLUSIONS 

In this work, we enumerated and classified the mechan- 
ically stable packings of bidisperse frictionless disks that 
occur as a function of the applied shear strain 7. We 
showed that MS packings form continuous geometrical 
families defined by the network of particle contacts. 

In addition, we studied the evolution of these systems 
during quasistatic shear strain at zero pressure to mimic 




y 

FIG. 22: Evolution of <f)j during quasistatic shear flow for 
N = 14 for a system with (dashed line) and without (solid 
line) Gaussian noise with width 0.01 times the small particle 
diameter added to the particle positions after each strain step. 

the dynamics of slowly sheared granular media. For small 
systems N < 16, we found that the dynamics was deter- 
ministic and strongly contracting, i.e. if the system is 
initialized in an ensemble of MS packings, it will quickly 
contract to at most a few MS packing families. The 
strong contraction stems from an abundance of short 
families, a propensity for the system to undergo more 
jumps than kinks in ^7(7), the fact that jumps only lead 
to increases in packing fraction, and the observation that 
families at higher packing fraction attract more jumps. 

In our studies of system sizes N > 16, we began to 
see features of large systems, including a dramatic in- 
crease in the period of the periodic orbits and bifurcations 
that lead to the random splitting of trajectories (</>,/ vs. 
7). We suggest that both the contraction and splitting 
mechanisms will persist in the large-system limit, and the 
fraction of MS-packing geometrical families that are vis- 
ited in steady-state will depend on ratio of the splitting 
and contraction rates. In large systems, we suspect that 
the dynamics will focus the system onto sets of frequent 
MS-packing families with similar structural and mechan- 
ical properties, although much more work is required to 
quantify these claims. 

Our long-term research objective is to develop a 
master-equation formalism to describe macroscopic 
slowly driven granular systems from the 'bottom-up' in 
terms of collections of small subsystems or microstates. 
In this manuscript, we took a significant step forward in 
this effort. We identified the types of microstates that 
can exist over the full range of shear strain and studied 
the probabilities with which they occur. This information 
can be used as input in the master-equation approach 
to calculate the contraction and splitting rates and ulti- 
mately the steady-state distributions of macroscale MS 
packings. 
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APPENDIX A: NUMERICAL DETAILS 



TABLE V: The eight roto-inversion matrices R; in 2D that 
when applied to a given MS packing at integer values of shear 
strain generate the eight equivalent polarizations with identi- 
cal eigenvalue lists of the dynamical matrix. 




In this appendix, we elaborate on technical details 
of the simulations not described in the main text. We 
include specific numerical parameters of the packing- 
generation protocol and the method to construct the 
complete geometrical family map. 



Packing- generation protocol In Sec. lIV Aj we outlined 
our procedure to generate mechanically stable packings. 
Here, we provide some of the numerical parameters in- 
volved in the simulations. For the energy minimization, 
we employ the conjugate gradient technique [32j, where 
the particles are treated as massless. The stopping crite- 
ria for the energy minimization (Vt — Vt-\ < Vtoi = 10~ 16 
and V t < Vmin = 10~ 16 , where V t is the potential 
energy per particle at iteration t) and the target po- 
tential energy per particle of a static granular packing 
(Vt.ni < V < 2Vtoi) are the same as used in previous stud- 
ies [17( • For the first compression or decompression step 
we use the packing-fraction increment Ac/) = 1CP 4 . Each 
time the procedure switches from expansion to contrac- 
tion or vice versa, A<j> is reduced by a factor of 2. Using 
the packing generation procedure with these parameters, 
we are able to locate the jamming threshold in packing 
fraction 4>j to within 10 -8 for each static packing over the 
full range of 7. Since we implement an energy minimiza- 
tion technique with no inertia, we do not need to alter 
the stopping criteria to handle rattler particles, which 
possess fewer than three contacts and are not members 
of the force bearing network. 

Geometrical Families To construct the complete map 
of geometrical families, we divided the region 7 = [0, 0.5] 
into small shear strain intervals 7 s +i — j s = A7 = 10~ 2 . 
For the range of system sizes N = 4 to 20 studied, this 
choice for A7 limited the number of rearrangement events 
to roughly one per shear strain interval. At each sampled 
shear strain j s , we generated at least N t = 10 6 MS pack- 
ings using random initial particle positions. 

Two MS packings at different shear strains are consid- 
ered to belong to the same geometrical family if they pos- 
sess the same set of particle contacts. The particle con- 
tact networks of two MS packings can be distinguished 
by comparing the eigenvalues of their connectivity ma- 
trices Cij , where the i j-th element of C is 1 if particles i 
and j are in contact and otherwise. Two systems have 
the same contact network if all of the eigenvalues of their 
connectivity matrices are the same. 



APPENDIX B: DYNAMICAL MATRIX 

In this appendix, we calculate the elements of the dy- 
namical matrix @ for the repulsive linear spring poten- 
tial ([8]). We employ slightly different notation for the dy- 
namical matrix compared to ^ by separating the spatial 
a, [3 = x, y and particle i,j = 1, ■ • ■ , N indexes. As given 
in Ref. [341 ] . for pairwise, central potentials the dynam- 
ical matrix has the following form for the off-diagonal 
components i j: 

Mia,j{3 = —— (S a - hjahjp) ~ CijhjahjP, (Bl) 
r ij 

where f« a is the ath component of r^ , 



dV{n 

dri. 



and 



_ d 2 V(r tJ ) 



cV 2 - 



9 



= - e _i 



(B2) 



(B3) 



In the calculation of tij and c^ , we have ignored the 6- 
function contributions arising from cases when particles 
i and j are just touching with r^ — o~ij. The diagonal 
components (i — j) of the dynamical matrix are given by 



N 



M, 



^2 M ia,]P- 



(B4) 



The shear-periodic boundary conditions only affect the 
definition of the separation for particles near the edges 
of the simulation cell. 



APPENDIX C: POLARIZATIONS OF MS 
PACKINGS 

In this appendix, we provide details for generating dif- 
ferent polarizations in simulations and determining which 
polarizations are distinct. 

We will first describe the symmetries that MS packings 
possess under shear periodic boundary conditions in un- 
deformed cells, since these symmetries affect the number 
and types of MS packings that can be obtained during 
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(a) (b) 




(c) (d) 




FIG. 23: (a) Typical MS packing at half-integer shear strains. 
We apply the transformations that are consistent with shear- 
periodic boundary conditions — the roto-inversion transforma- 
tions R3, Re, and Rg in Table IVl — to the configuration in (a) 
to generate the configurations in (b)-(d). However, (a) and 
(b) and (c) and (d) are related by the shear symmetry opera- 
tion (rotations by it, R3), and thus there are only two distinct 
polarizations at half integer strains. 

shear. At integer shear strains, there are eight 'polar- 
izations' all of which have the same list of eigenvalues 
of the dynamical matrix, but different eigenvectors. A 
given MS packing at integer shear strain (panel (1)) and 
its equivalent polarizations (panels (2)-(8)) are shown in 
Fig. O Each polarization with coordinates £j in panel (i) 
in Fig. [5] can be obtained from the original MS packing 
£1 in panel (1) using 

6 = Rili, (Ci) 



where the eight roto-inversion matrices Ri with detRj = 
±1 in 2D are given in Table [Vj In isotropically com- 
pressed systems these polarizations occur with equal 
probability. However, in systems subjected to shear 
strain, these polarizations will occur with different prob- 
abilities. 

2D systems subjected to simple planar shear flow pos- 
sess a discrete rotational symmetry; i.e. the system is 
unchanged when it is rotated by it about an axis coming 
out of the page (i.e. apply R3 to a given MS packing) 
as shown in the bottom of Fig. [5] In panels (l)-(8) in 
Fig. [5j we see that polarizations 1 and 5, 2 and 6, 3 and 
7, and 4 and 8 are related by rotations by 7t, and therefore 
will behave the same under simple planar shear. Thus, 
only four distinct polarizations at integer shear strains 
remain. 

If the accumulated shear strain is half-integer, MS 
packings will have at most only two distinct polariza- 
tions. Only four roto-inversion transformations are con- 
sistent with shear-periodic boundary conditions, R3, R6, 
and Rg in Table |V1 These transformations have been ap- 
plied to the configuration (a) in Fig. [23] to generate the 
configurations (b)-(d). However, (a) and (b) and (c) and 
(d) are related by the shear symmetry operation (rota- 
tions by 7r, R3), and thus only two distinct polarizations 
remain at half integer strains. 

In general, there is only one polarization for all shear 
strains other than integer and half-integer. Moreover, the 
number of polarizations at integer (half-integer) strains 
can be smaller than four (two) if particular MS packings 
possess additional symmetries. 

In simulations, we distinguish the polarizations of 
two MS packings by applying all possible roto-inversion 
transformations configurations consistent with the shear- 
periodic boundary conditions to a given MS packing. 
We then compare the eigenvector (corresponding to the 
smallest nonzero eigenvalue) of the second configuration 
to that for all of the roto-inverted configurations of the 
first and look for a match. 
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